Fourier's law from a chain of coupled anharmonic oscillators under energy conserving 
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We analyze the transport of heat along a chain of particles interacting through anharmonic po- 
tentials consisting of quartic terms in addition to harmonic quadratic terms and subject to heat 
reservoirs at its ends. Each particle is also subject to an impulsive shot noise with exponentially 
distributed waiting times whose effect is to change the sign of its velocity, thus conserving the en- 
ergy of the chain. We show that the introduction of this energy conserving stochastic noise leads 
to Fourier's law. The behavior of thels heat conductivity for small intensities of the shot noise 
and large system sizes are found to obey a finite-size scaling relation. We also show that the heat 
conductivity is not constant but is an increasing monotonic function of temperature. 
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I. INTRODUCTION 

The derivation of Fourier's law, or any other macro- 
scopic law, from the microscopic underlying motion of 
particles constitutes a major task in condensed matter 
physics. This task comprises not only the derivation it- 
self but also the problem of setting up the appropriate 
microscopic model. The simplest model one could con- 
ceive to derive Fourier's law is a chain of particles in- 
teracting through harmonic potentials in contact with 
two heat reservoirs at each end. However, it has been 
shown by Rieder et al [l| that this model does not lead 
to Fourier's law. Since then, several microscopic models 
have been introduced and studied (§~22] some of them 
leading instead to the so called anomalous Fourier's law. 

Fourier's law states that the heat flux J is proportional 
to the gradient of the temperature T, that is, 



dT 

J = -K-, 
dx 



(1) 



where n is the heat conductivity. If we consider a small 
bar of length L subject to a difference in temperature 
AT, then J = kAT/L. Thus a microscopic model for 
Fourier's law should predict a heat flux that decreases 
with L, for a fixed value of AT, according to 
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(2) 



This amounts to saying that k is finite. If, instead, we 
find that J ~ L~ a with a < 1, then we are faced with 
the anomalous Fourier's Law. In this case we may say 
that k is infinite, as in the harmonic chain diverging 
according to k ~ L a , with exponent a = 1 — a. Notice 
that L should be macroscopically small, so that equation 
@ is the expression of Fourier's law, but microscopically 
large so that a microscopic model for this law should yield 
equation ^ for sufficiently large L. 

The heat flux J of the linear harmonic chain placed 
between two heat reservoirs has been shown [l[ to be in- 
dependent of the size L of the chain, meaning that the 
heat conductivity k is infinite. This result is a direct 
consequence of the ballistic transmission of heat by the 



elastic waves, from one reservoir to the other. A conse- 
quence of this result is that a perfectly harmonic crystal 
has an infinite heat conductivity [24j . In real crystals the 
heat conductivity is manifestly finite. This is due to the 
presence of lattice imperfections, impurities, and other 
factors that act as scattering centers for the waves carry- 
ing energy, such as the Umklapp process [25[ . These fac- 
tors make the heat conduction a diffusion process which 
implies Fourier's law. The crossover from ballistic do dif- 
fusive behavior can also be undertood in terms of the 
dephasing of the elastic waves caused by the scattering 
events just mentioned. This process is characterized by 
a dephasing length [2(| [53] which, when smaller than the 
size of the system, gives the linear temperature profile ex- 
pected from Fourier's law. Thus, a possible ingredient in 
the microscopic derivation of Fourier's law consists in the 
presence of a diffusive motion at the microscopic level. 

One way of introducing this ingredient is by means of 
stochastic collisions that change the sign of the velocity 
of the particles. This can be accomplished in the form 
of impulsive shot noises with exponentially distributed 
waiting times (22j . Two key properties are required when 
devising such a noise. Firstly, it should conserve the total 
energy of the system because any variation of the energy 
of the system should only be due to the contact with 
the heat reservoirs. Changing the sign of the velocity 
does not alter the energy. Secondly, it should make the 
system ergodic even when it is not coupled to the heat 
baths. A harmonic chain with this type of shot noise has 
been indeed studied by Dhar et al 22] . They showed that 
this model can be mapped into the self consistent har- 
monic crystal introduced by Bolsteri et al [28| , and stud- 
ied by Bonetto et al [ll|- In this model each particle is 
in contact with independent heat reservoirs, whose tem- 
peratures are chosen so that, in the steady state, there 
is no exchange of heat between these reservoirs and the 
chain. The contact with the reservoirs is regarded as a 
procedure to make the system ergodic (28|. This model 
predicts a linear profile for the temperature and a finite 
heat conductivity, which is independent of temperature. 

Here we study a chain of particles interacting through 
anharmonic forces in addition to random reversals of 
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the velocity. More specifically we consider a potential 
with quartic terms in addition to the harmonic quadratic 
terms. Without the stochastic shot noise, this is the 
well known Fermi-Pasta-Ulam model in contact with heat 
reservoirs at its ends, which was studied numerically by 
Lepri et al. Q who found a superdiffusive transport of 
heat implying an anomalous Fourier's law. This impor- 
tant result has been confirmed by other numerical stud- 
ies and other approaches 0, i, Qj, M, M, M, M, HI HI 
The impulsive stochastic shot noise that we use here can 
be regarded as a procedure that turns the superdiffusive 
transport of heat into a diffuse transport leading thus to 
Fourier's law. 

By numerically solving the Langevin equations for 
chains of several sizes, we determine the heat conduc- 
tivity as a function of the system size L and the rate of 
stochastic collisions A. When A is nonzero, we obtain a 
finite heat conductivity and therefore Fourier's law. For 
small values of A, the heat conductivity is found to be- 
have as as k ~ A~ b , diverging when A = 0. Our numeri- 
cal results give 6 = 0.52±0.06. We have also determined 
the exponent a related to the divergence of n with L at 
A = and found a = 0.42 ± 0.04. These results arc 
distinct from the harmonic case [l|, [22J for which a = 1 
and b = 1. Also, in contrast to the harmonic case, we 
have found that the heat conductivity depends on tem- 
perature. More precisely, for a fixed AT, we found that 
it increases monotonically with the temperatures of the 
reservoirs. 



II. MODEL 
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FIG. 1: Average squared velocity (v^) of the particles as a 
function of the position n on the chain, which is in contact 
with heat reservoirs at temperatures Ta = 2 and Tb = 1. 
Results are given for a chain of L = 40 particles for the values 
of A indicated, for the (a) harmonic case and (b) anharmonic 
case. In (a) the continuous lines are exact results. 



We consider a chain of L interacting particles with 
equal masses m and denote by x n the position of the 
n-th particle. The total potential energy V(x) = 
V(xi, X2, • • • , Xl) is considered to be a sum of anhar- 
monic potential energies involving nearest-neighbor pairs 
of particles 



V(x) = 



n=0 



x n+ if + -p(a;. 



(3) 

where K\ and Ky. are parameters. We consider fixed 
boundary conditions so that xl+i — xo = 0. When 
K 2 = 0, the harmonic potential is recast. The force F n 
acting on the n-th particle due to the potential V(x) is 

F n (x) = Ki(x n -i + x n +i - 2x n )+ 



K 2 {x r , 



i) 3 + K2{Xn + l - Xnf 



(4) 



The first and last particles are connected to heat baths 
at temperatures Ta and Tb, and all particles are sus- 
ceptive to stochastic collisions, here described by forces 
T n (t). Denoting by v n = dx n /dt the velocity of the n-th 
particle, the equations of motion are stochastic equations 



given by 



m lt =Fl+ ~ aVl + ^ 2ak B T AU, (5) 



dv 

m-£- = F n + F n {t), 2 < n < L - 1, (6) 



m ~dT = Fl+ W _ aVL + V 2a kBT B ^B, (7) 

where fcs is Boltzmann constant and the last two terms 
in equations (0 and (0 represent the coupling to the 
heat reservoirs. The constant a represents the strength 
of the coupling and £,a and £b are independent Gaus- 
sian white noises with zero mean and unit variance. The 
forces !F n {t) have the form of impulsive shot noises, given 
by 



J n 



(8) 



where t n g, £ — 1,2,... are uncorrelated exponentially dis- 
tributed stochastic waiting times with a probability den- 
sity distribution p{t) — Ae~ At . Here, the parameter A 
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is the rate of collisions, which has been taken to be the 
same for all particles. After a collision occurring at time 
t, the n-th particle changes its velocity from v n (t~) to 
v n (t + ) = —v n (t~), thus conserving its kinetic energy and 
therefore the total energy 

E = jr^vl + V{x). (9) 

71=1 

Due to the contact with the reservoirs, the total en- 
ergy E is not a strictly conserved quantity. From the 
stochastic equations of motion we find 

^- = Ja + Jb, (10) 
at 

where 

Ja = -av\ + v 1 y / 2ak B T A £,A, (11) 

and 

J B = -av\ + v L ^2ak B T B £ B . (12) 

In the stationary state, (E) is a constant and the sum of 
the fluxes Ja = {Ja) and J B = (Jb) vanishes, that is, 
J B = — J a- The heat flux J = J a can be calculated as 
the average of the right-hand side of equation (ITT]) or in 
a equivalent way from the equation 

J = Ki((x n -i - x n )v n ) + K 2 ((x n -i - x n ) 3 v n ), (13) 

which we found to be numerically more accurate than 
the formula J = {Ja)- 

III. NUMERICAL SOLUTIONS 
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FIG. 2: Heat flux J as a function of the inverse of the system 
size L for several values of A for (a) the harmonic chain and (b) 
the anharmonic chain. In both cases, from top to bottom, A = 
0.05,0.1,0.2,0.5 and 1. The temperatures of the reservoirs 
are Ta = 2 and Tb = 1. In (a) the continuous lines are exact 
results. 



The stochastic equations of motion were solved nu- 
merically for chains of several sizes L. This was done 
by discretizing the time in intervals At. We use an ap- 
proach in which the deterministic part of the equations 
of motion of the inner particles are handled by the use 
of the Verlet algorithm [30] so as to ensure that, in the 
absence of the heat baths, energy is conserved. For the 
equations of motion of the first and last particles, which 
contain the stochastic forces due to the heat baths, we 
used the stochastic Verlet algorithm developed in |3l| . 
As for the stochastic shot noises we treat them as fol- 
lows. At each time step, the sign of the velocity of each 
particle is changed with probability p = XAt. This pro- 
cedure generates a Poisson process with discrete waiting 
time t = lAt that is distributed according to the proba- 
bility distribution p(l — p) e . In the continuous time limit, 
At 0, this yields the exponential distribution Ae _At , 
as required. 

For definiteness, our numerical calculations were per- 
formed with k B — 1, m — 1, a — 1 and At = 0.01. For 
the anharmonic potential all results reported in this pa- 
per were obtained for K\ — 1 and K 2 = 1. The size of 



the system ranged from L = 10 up to L = 5000. We also 
present numerical results for the harmonic case (K 2 = 0), 
for K\ = 1 and compare with the results of the anhar- 
monic case. The existing exact solution [22j for harmonic 
case is used to check our numerical procedure. The ex- 
act solution is obtained by solving the equations for the 
pair correlations, which is possible because they consist 
of closed equations. However, this closure property does 
not happen for the anharmonic case. 

In Fig. [JJ we show the results for the average kinetic 
energy for each particle as a function of the position n on 
the chain for the harmonic and anharmonic cases. The 
temperatures of the reservoirs are considered to be dis- 
tinct, Ta ^ T B , and our numerical calculations were per- 
formed for Ta = 2 and T B = 1 . Without stochastic colli- 
sions (A = 0) the results of the harmonic case shows that 
the kinetic energy is almost constant, a result obtained 
by Rieder et al [Ij, which does not lead to Fourier's Law. 
The inclusion of the stochastic collisions (A / 0) pro- 
duces a drastically different result. The average kinetic 
energy as a function of n displays now a nonzero slope 
as can be seen in Fig. [TJ For the anharmonic case, all 
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FIG. 3: Log-log plot of the heat conductivity Hi ELS Su function 
of the system size L for (a) the harmonic chain and (b) the an- 
harmonic chain. From top to bottom, A = 0, 0.02, 0.05, 0.1, 0.2 
and 0.5, for (a) and A = 0,0.001,0.002,0.005,0.01,0.02,0.05 
and 0.1, for (b). The temperatures of the reservoirs are 
Ta = 2 and Tb = 1. In (a), the continuous lines are ex- 
act results and the slope of the straight line corresponding to 
A = is a = 1. In (b), the slope of the straight line (dashed 
line) fitted to the data points with large L corresponding to 
A = gives a = 0.42. 
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FIG. 4: Log-log plot of the reciprocal of the heat conductivity 
function of the collision noise A for (a) the harmonic 
chain and (b) the anharmonic chain, for several values of L 
indicated. The temperatures of the reservoirs are Ta = 2 
and Tb = 1. In (a) the continuous lines are exact results. 
The slope of the straight dashed lines is b = 1 for (a) the 
harmonic case and b = 0.52 for the (b) anharmonic case. The 
lowest curve in both cases are extrapolation obtained from 
finite values of L. 



curves, including the case A = 0, show a nonzero slope. 
In spite of that, the A = case does not lead to Fourier's 
law. 

We have calculated the flux J at the stationary state 
by using equation (|T5)> and the results are shown in Fig. 
H as a function of 1/L for several values of A. From 
this figure we see clearly that J ~ 1/L, for sufficiently 
large values of L, in accordance with Fourier's law, as 
long as A ^ 0, for both harmonic and anharmonic cases. 
For sufficiently large values of L the heat conductivity 
is given by n — JL/AT. This quantity is plotted as a 
function of L for several values of the rate of stochastic 
collisions A, including A = 0. The results are shown in 
Fig. [3] For both the harmonic and anharmonic cases, 
the heat conductivity k approaches a constant when L — > 
oo, as long as A ^ 0, and Fourier's law is accomplished. 
When A = 0, our numerical results gives a superdiffusive 
behavior with k —> oo when L — > oo, according to 

k ~ L a , A = 0, (14) 

as shown in Fig. [3] For the harmonic case, a = 1, which 
is in accordance with the result by Rieder et al [l|. For 



the anharmonic case we get a = 0.42 ± 0.04. This is in 
agreement with the result a = 0.45 ± 5 found by Lepri et 
al. J2| and in excelent agreement with the value a = 2/5 

BQ3. 

To analyze the behavior of the heat conductivity k = 
JL/AT as A 4 we have plotted this quantity as a 
function of A for several values of the system size L. The 
results are shown in Fig. |4j We have plotted also the 
extrapolated values of k when L —> 00 for each A. These 
extrapolated values were extracted from the plot of 1/k 
versus l/L. The heat conductivity k diverges when A — > 
according to 

k ~ X~ b , L^oo. (15) 

For the harmonic case our results give b = 1, a result 
obtained by Dhar et al [22[ • For the anharmonic case we 
found b = 0.52 ± 0.06, a result clearly distinct from the 
harmonic case. 

The algebraic behavior of n with L and A can be ob- 
tained by assuming the following finite-size scaling for 
the heat conductivity, 

k = L a vjj(\L c ), (16) 
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where ip(s) is a universal function of s = XL C such that 
ip(Q) is a finite constant and ip(s) ~ s~ b when s is large. 
To ensure a finite conductivity in the limt L — > oo, the 
exponent c must be related to the exponents a and 6 by 
c = a/6. In Fig. [5j we show the data collapse for the 
harmonic case by plotting n/L a as a function of XL C , 
where in this case a = 1 and c = 1. The data collapse is 
well described by the expression 

M = fAr, (it) 

as can be seen in Fig. [5] This function was obtained 
by numerically solving the exact equations for the pair 
correlations from which we found A = 1/4 and B = 3/2. 
For large values of s, this result gives -0 ~ s _1 so that 
k ~ A -1 , independent of L, which is the behavior of 
the conductivity for the harmonic chain when L — > oo 
[llT [22l [28| . For the anharmonic case, the data colapse, 
shown in Fig. [5j was obtained by using the exponents a 
and 6 obtained previously. 

It is worth mentioning that Lepri et al. [32| in their 
study of a one-dimensional harmonic crystal with energy 
conservative noise consisting of elastic collisions between 
neighboring particles reported a conductivity that be- 
haves as (L/A) 1//2 for large values of L. Assuming for 
this system a scaling function of the form (fTB| and keep- 
ing in mind that a = 1 , because the conductivity behaves 
as k ~ L when A = 0, we may conclude that 6 = 1/2 and 
c = 1. Notice, however, that c ^ a/6 so that n is not 
finite but diverges when L — Y 00. 

A relevant feature of the present anharmonic chain 
with random reversal of velocities is the dependence of 
the heat conductivity with temperature. For the har- 
monic case the heat conductivity is temperature indepen- 
dent [22], [28| , a result that can be understood by using 
the following reasoning. If we rescale the temperature 
of the reservoirs by a factor r, that is, T A — > rT A and 
Tb — > tTb, and the positions and velocities by a factor 
r 1 / 2 , that is, v„ — > r x / 2 w n and x n — > r l l 2 x n , the equa- 
tions of motions for the case Ki = become invariant. 
In addition, according to Eq. (jT3"]) with K2 = 0, the heat 
flux changes as the temperature, that is, J — > rJ. From 
this relation it follows that J is a homogeneous func- 
tion of T A and T B , so that J = T A cj)(T B /T A ). Writing 
T B = T A + AT, then for small values of AT it follows 
that J ~ AT leading to a heat conductivity k indepen- 
dent of temperature, a result that we have also checked 
numerically. 

For the anharmonic case we have found that the heat 
conductivity k is an increasing function of temperature. 
In Fig. |6j we show k as a function of the temperature 
T = T B of the colder reservoir, The heat conductivity 
was determined from n — JL/ AT, for the same value of 
AT = 1 and for several values of the noise parameter A. 
We used L = 500, a value high enough so that the values 
of k may be considered the asympototic values (see Fig. 
13]), with the exception of the case A = 0. Our results 
indicate a monotonic increase of the heat conductivity 
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FIG. 5: Data colapse obtained from the numerical results of 
the heat conductivity k obtained from several values of A and 
L corresponding to the harmonic case (a) and the anharmonic 
case (b), where = njL a and s — \L a ^ b . The continuous line 
in case (a) is described by the expression (|17|) with A = 1/4 
and B — 3/2. In case (b), the continous line is a guide to 
the eyes. The temperatures of the reservoirs are Ta = 2 and 
Tb ~ 1. Points with the same symbol correspond to the same 
value of A. 

with temperature as can be seen in Fig. [S] We have 
found that our results are consistent with the result of 
Aoki and Kusnezov 1331 and with the upper bounds of 
Bernardin and Olla [34j |. 



IV. CONCLUSION 

In conclusion, we have considered a chain of particles 
interacting through anharmonic potentials and subject to 
heat reservoirs at its ends. In addition, the chain is sub- 
ject to a shot noise that changes the sign of the velocities 
of the particles at random times, distributed according 
to an exponential distribution. The shot noise does not 
change the energy so that the changes in energy are only 
due to the contact with the reservoirs. We have shown 
that, when the chain is connected to reservoirs at differ- 
ent temperatures, the heat flux is inversely proportional 
to the size of the system, as long as the shot noise pa- 
rameter A is nonzero, and therefore in accordance with 
Fourier's law. Our results suggest, in accordance with 
[28|, that the ergodicity may play a crucial role in the 
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FIG. 6: Heat conductivity k = JL/AT as a function of the 
temperature T = Tb of the colder reservoir, for several values 
of A. From top to bottom A = 0, 0.02, 0.05 and 0.1. All results 
were obtained for a chain of size L = 500 and by using the 
same value of the difference in temperature of the reservoirs, 
AT — 1. The continuous lines are guide to the eyes. 



derivation of Fourier's law. 

We have also obtained the behavior of k with L when 
A = and k with A when L — » oo. Both behaviors 
are found to be algebraic, characterized by exponents 
a = 0.42 ± 0.04 and b = 0.52 ± 0.06. This allows us 
to introduce a finite-size scaling in which the noise pa- 
rameter scales with the inverse of the system size to the 
power a = b/c. For the harmonic case, we have shown 
that a finite size scaling exists, a result that has been 
also obtained by numerically solving the equations for 
the correlations. 
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